The question is how much G x E we see in the RNAseq data.

To get started on this I ran the file “ProcessCount.R”. For this the data was TMM normalized, voom transformed, and then a gtXenvironment model was fit in limma.

library(limma)
library(tidyverse)
── Attaching packages ────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 2.2.1     ✔ purrr   0.2.4
✔ tibble  1.3.4     ✔ dplyr   0.7.4
✔ tidyr   0.7.2     ✔ stringr 1.2.0
✔ readr   1.1.1     ✔ forcats 0.2.0
── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(stringr)
library(modelr)
library(broom)

Attaching package: ‘broom’

The following object is masked from ‘package:modelr’:

    bootstrap
library(multidplyr)
library(VCA)

Get the limma objects

load("voom.int.Rdata")

How many genes show GxE?

generate vector of interaction coefficients and then topTable

colnames(fit.int$coefficients)
int.coef.n <- colnames(fit.int$coefficients) %>% str_which(":") 
topTable(fit.int,coef = int.coef.n)

Not working

let’s just do an ANOVA for each gene

counts <- as.tibble(counts.voom.gr$E) %>% 
  mutate(geneIndex=1:n())
counts[1:10,1:10]
counts.long <- counts %>% gather(key="name",value="expression", -geneIndex) %>%
  mutate(gt=factor(str_extract(name,"RIL_[0-9]*")), 
         env=factor(str_extract(name,"CR|UN")), 
         rep = factor(str_extract(name,"Rep[0-9]"))) %>%
  select(-name) %>%
  group_by(geneIndex) %>% nest()
counts.long
counts.long$data[[1]] 
#test.df <- counts.long$data[[1]] 

Now fit an anova and do a VCA for each gene in parallel

fitAnova <- function(df) aov(expression ~ gt*env,data=df)

getPvals <- function(fit) {
  pvals <- summary(fit)[[1]]$`Pr(>F)`
  names(pvals) <- dimnames(summary(fit)[[1]])[[1]] %>% str_trim()
  data.frame(t(pvals))
}

fitVCA <- function(df) {
  tmpvca <- anovaVCA(expression ~ gt*env,Data=as.data.frame(df))
  data.frame(t(tmpvca$aov.tab[,"%Total"]))
}

cl <- create_cluster(cores=4)
counts.long <- partition(counts.long,geneIndex,cluster = cl)
cluster_library(cl,c("purrr","broom","VCA","stringr"))
cluster_copy(cl,fitAnova)
cluster_copy(cl,getPvals)
cluster_copy(cl,fitVCA)

system.time(
  counts.results <- counts.long %>% mutate(aov = map(data,fitAnova), 
                                        glance = map(aov,glance),
                                        pvals = map(aov,getPvals),
                                        varcomp = map(data,fitVCA)) %>%
    select(-data, -aov)
)
system.time(cluster_rm(cl,"counts.long"))
rm(counts.long)
system.time(counts.results <- collect(counts.results))
save(counts.results,file="GxE_models.Rdata")
(aovsummary <- counts.results %>% unnest(glance,.drop = TRUE))
(pvals <- counts.results %>% unnest(pvals,.drop = TRUE))
(varcomp <- counts.results %>% unnest(varcomp,.drop = TRUE))

PVals

Number of genes (out of 24619) that have a FDR adjusted pvalue of less than 0.05:

pvals %>% ungroup() %>%
  summarise(gt=sum(p.adjust(gt,method="fdr") < 0.05),
            env=sum(p.adjust(env) < 0.05),
            gt.env=sum(p.adjust(gt.env,method="fdr") < 0.05))

Variance components

Histograms of % variance explained

pl <- varcomp %>%
  ungroup() %>%
  select(-total) %>%
  gather(key = "component", value = "percent.variance",-geneIndex) %>%
  ggplot(aes(x=percent.variance)) +
  facet_wrap(~ component, scale="free_y") +
  geom_histogram()
pl

focus in on gt x env histogram, non-zeros

pl <- varcomp %>%
  ungroup %>%
  filter(gt.env > 0) %>%
  ggplot(aes(x=gt.env)) +
  geom_histogram(binwidth = 5)
pl

Number of genes with gt.env explaining X % variance:

varcomp$gt.env %>%
  cut(breaks=seq(0,35,5),right = FALSE) %>% table
.
  [0,5)  [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) 
  23311    1121     151      29       5       1       1 
LS0tCnRpdGxlOiAiRyB4IEUgZXhwcmVzc2lvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhlIHF1ZXN0aW9uIGlzIGhvdyBtdWNoIEcgeCBFIHdlIHNlZSBpbiB0aGUgUk5Bc2VxIGRhdGEuICAKClRvIGdldCBzdGFydGVkIG9uIHRoaXMgSSByYW4gdGhlIGZpbGUgIlByb2Nlc3NDb3VudC5SIi4gIEZvciB0aGlzIHRoZSBkYXRhIHdhcyBUTU0gbm9ybWFsaXplZCwgdm9vbSB0cmFuc2Zvcm1lZCwgYW5kIHRoZW4gYSBndFhlbnZpcm9ubWVudCBtb2RlbCB3YXMgZml0IGluIGxpbW1hLiAgCgpgYGB7cn0KbGlicmFyeShsaW1tYSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc3RyaW5ncikKbGlicmFyeShtb2RlbHIpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkobXVsdGlkcGx5cikKbGlicmFyeShWQ0EpCmBgYAoKR2V0IHRoZSBsaW1tYSBvYmplY3RzCmBgYHtyfQpsb2FkKCJ2b29tLmludC5SZGF0YSIpCmBgYAoKSG93IG1hbnkgZ2VuZXMgc2hvdyBHeEU/CgpnZW5lcmF0ZSB2ZWN0b3Igb2YgaW50ZXJhY3Rpb24gY29lZmZpY2llbnRzIGFuZCB0aGVuIHRvcFRhYmxlCmBgYHtyLCBldmFsPUZBTFNFfQpjb2xuYW1lcyhmaXQuaW50JGNvZWZmaWNpZW50cykKaW50LmNvZWYubiA8LSBjb2xuYW1lcyhmaXQuaW50JGNvZWZmaWNpZW50cykgJT4lIHN0cl93aGljaCgiOiIpIAp0b3BUYWJsZShmaXQuaW50LGNvZWYgPSBpbnQuY29lZi5uKQpgYGAKCk5vdCB3b3JraW5nCgpsZXQncyBqdXN0IGRvIGFuIEFOT1ZBIGZvciBlYWNoIGdlbmUKCmBgYHtyfQpjb3VudHMgPC0gYXMudGliYmxlKGNvdW50cy52b29tLmdyJEUpICU+JSAKICBtdXRhdGUoZ2VuZUluZGV4PTE6bigpKQpjb3VudHNbMToxMCwxOjEwXQpgYGAKCmBgYHtyfQpjb3VudHMubG9uZyA8LSBjb3VudHMgJT4lIGdhdGhlcihrZXk9Im5hbWUiLHZhbHVlPSJleHByZXNzaW9uIiwgLWdlbmVJbmRleCkgJT4lCiAgbXV0YXRlKGd0PWZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJSSUxfWzAtOV0qIikpLCAKICAgICAgICAgZW52PWZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJDUnxVTiIpKSwgCiAgICAgICAgIHJlcCA9IGZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJSZXBbMC05XSIpKSkgJT4lCiAgc2VsZWN0KC1uYW1lKSAlPiUKICBncm91cF9ieShnZW5lSW5kZXgpICU+JSBuZXN0KCkKY291bnRzLmxvbmcKY291bnRzLmxvbmckZGF0YVtbMV1dIAojdGVzdC5kZiA8LSBjb3VudHMubG9uZyRkYXRhW1sxXV0gCmBgYAoKTm93IGZpdCBhbiBhbm92YSBhbmQgZG8gYSBWQ0EgZm9yIGVhY2ggZ2VuZSBpbiBwYXJhbGxlbApgYGB7ciwgZXZhbD1GQUxTRX0KZml0QW5vdmEgPC0gZnVuY3Rpb24oZGYpIGFvdihleHByZXNzaW9uIH4gZ3QqZW52LGRhdGE9ZGYpCgpnZXRQdmFscyA8LSBmdW5jdGlvbihmaXQpIHsKICBwdmFscyA8LSBzdW1tYXJ5KGZpdClbWzFdXSRgUHIoPkYpYAogIG5hbWVzKHB2YWxzKSA8LSBkaW1uYW1lcyhzdW1tYXJ5KGZpdClbWzFdXSlbWzFdXSAlPiUgc3RyX3RyaW0oKQogIGRhdGEuZnJhbWUodChwdmFscykpCn0KCmZpdFZDQSA8LSBmdW5jdGlvbihkZikgewogIHRtcHZjYSA8LSBhbm92YVZDQShleHByZXNzaW9uIH4gZ3QqZW52LERhdGE9YXMuZGF0YS5mcmFtZShkZikpCiAgZGF0YS5mcmFtZSh0KHRtcHZjYSRhb3YudGFiWywiJVRvdGFsIl0pKQp9CgpjbCA8LSBjcmVhdGVfY2x1c3Rlcihjb3Jlcz00KQpjb3VudHMubG9uZyA8LSBwYXJ0aXRpb24oY291bnRzLmxvbmcsZ2VuZUluZGV4LGNsdXN0ZXIgPSBjbCkKY2x1c3Rlcl9saWJyYXJ5KGNsLGMoInB1cnJyIiwiYnJvb20iLCJWQ0EiLCJzdHJpbmdyIikpCmNsdXN0ZXJfY29weShjbCxmaXRBbm92YSkKY2x1c3Rlcl9jb3B5KGNsLGdldFB2YWxzKQpjbHVzdGVyX2NvcHkoY2wsZml0VkNBKQoKc3lzdGVtLnRpbWUoCiAgY291bnRzLnJlc3VsdHMgPC0gY291bnRzLmxvbmcgJT4lIG11dGF0ZShhb3YgPSBtYXAoZGF0YSxmaXRBbm92YSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZ2xhbmNlID0gbWFwKGFvdixnbGFuY2UpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHZhbHMgPSBtYXAoYW92LGdldFB2YWxzKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHZhcmNvbXAgPSBtYXAoZGF0YSxmaXRWQ0EpKSAlPiUKICAgIHNlbGVjdCgtZGF0YSwgLWFvdikKKQpzeXN0ZW0udGltZShjbHVzdGVyX3JtKGNsLCJjb3VudHMubG9uZyIpKQpybShjb3VudHMubG9uZykKCmBgYAoKYGBge3IsIGV2YWw9RkFMU0V9CnN5c3RlbS50aW1lKGNvdW50cy5yZXN1bHRzIDwtIGNvbGxlY3QoY291bnRzLnJlc3VsdHMpKQpzYXZlKGNvdW50cy5yZXN1bHRzLGZpbGU9Ikd4RV9tb2RlbHMuUmRhdGEiKQpgYGAKCmBgYHtyfQppZiAoISBleGlzdHMoY291bnRzLnJlc3VsdHMpKSBsb2FkKCJHeEVfbW9kZWxzLlJkYXRhIikKKGFvdnN1bW1hcnkgPC0gY291bnRzLnJlc3VsdHMgJT4lIHVubmVzdChnbGFuY2UsLmRyb3AgPSBUUlVFKSkKKHB2YWxzIDwtIGNvdW50cy5yZXN1bHRzICU+JSB1bm5lc3QocHZhbHMsLmRyb3AgPSBUUlVFKSkKKHZhcmNvbXAgPC0gY291bnRzLnJlc3VsdHMgJT4lIHVubmVzdCh2YXJjb21wLC5kcm9wID0gVFJVRSkpCmBgYAoKIyMgUFZhbHMKTnVtYmVyIG9mIGdlbmVzIChvdXQgb2YgMjQ2MTkpIHRoYXQgaGF2ZSBhIEZEUiBhZGp1c3RlZCBwdmFsdWUgb2YgbGVzcyB0aGFuIDAuMDU6CmBgYHtyfQpwdmFscyAlPiUgdW5ncm91cCgpICU+JQogIHN1bW1hcmlzZShndD1zdW0ocC5hZGp1c3QoZ3QsbWV0aG9kPSJmZHIiKSA8IDAuMDUpLAogICAgICAgICAgICBlbnY9c3VtKHAuYWRqdXN0KGVudikgPCAwLjA1KSwKICAgICAgICAgICAgZ3QuZW52PXN1bShwLmFkanVzdChndC5lbnYsbWV0aG9kPSJmZHIiKSA8IDAuMDUpKQpgYGAKCiMjIFZhcmlhbmNlIGNvbXBvbmVudHMKCiMjIyBIaXN0b2dyYW1zIG9mICUgdmFyaWFuY2UgZXhwbGFpbmVkCgpgYGB7cn0KcGwgPC0gdmFyY29tcCAlPiUKICB1bmdyb3VwKCkgJT4lCiAgc2VsZWN0KC10b3RhbCkgJT4lCiAgZ2F0aGVyKGtleSA9ICJjb21wb25lbnQiLCB2YWx1ZSA9ICJwZXJjZW50LnZhcmlhbmNlIiwtZ2VuZUluZGV4KSAlPiUKICBnZ3Bsb3QoYWVzKHg9cGVyY2VudC52YXJpYW5jZSkpICsKICBmYWNldF93cmFwKH4gY29tcG9uZW50LCBzY2FsZT0iZnJlZV95IikgKwogIGdlb21faGlzdG9ncmFtKCkKcGwKYGBgCgpmb2N1cyBpbiBvbiBndCB4IGVudiBoaXN0b2dyYW0sIG5vbi16ZXJvcwpgYGB7cn0KcGwgPC0gdmFyY29tcCAlPiUKICB1bmdyb3VwICU+JQogIGZpbHRlcihndC5lbnYgPiAwKSAlPiUKICBnZ3Bsb3QoYWVzKHg9Z3QuZW52KSkgKwogIGdlb21faGlzdG9ncmFtKGJpbndpZHRoID0gNSkKcGwKYGBgCgpOdW1iZXIgb2YgZ2VuZXMgd2l0aCBndC5lbnYgZXhwbGFpbmluZyBYICUgdmFyaWFuY2U6CmBgYHtyfQp2YXJjb21wJGd0LmVudiAlPiUKICBjdXQoYnJlYWtzPXNlcSgwLDM1LDUpLHJpZ2h0ID0gRkFMU0UpICU+JSB0YWJsZQpgYGAKCg==